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Abstract. Landau-Zener-Stiickelberg (LZS) interference of continuously driven 
superconducting qubits is studied. Going beyond the second order perturbation 
expansion, we find a time dependent stationary population evolution as well 
as unsymmetrical microwave driven Landau-Zener transitions, resulting from the 
nonresonant terms which are neglected in rotating-wave approximation. For the 
low-frequency driving, the qubit population at equilibrium is a periodical function 
of time, owing to the contribution of the nonresonant terms. In order to obtain 
the average population, it is found that the average approximation based on 
the perturbation approach can be applied to the low-frequency region. For the 
extremely low frequency which is much smaller than the decoherence rate, we develop 
noncohcrcncc approximation by dividing the evolution into discrete time steps during 
which the coherence is lost totally. These approximations present comprehensive 
analytical descriptions of LZS interference in most of parameter space of frequency and 
decoherence rate, agreeing well with those of the numerical simulations and providing 
a simple but integrated understanding to system dynamics. The application of our 
models to microwave cooling can obtain the minimal frequency to realize effective 
microwave cooling. 
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1. Introduction 

Superconducting devises based on Joseplison tunnel junctions have great importance 
in demonstrating quantum phenomena at macroscopic scales and hold promise for 
applications in quantum information. [H El El HJ [5] In the studies of these 
systems, crossovers between energy levels are usually associated with a wide variety 
of interesting phenomena, particularly in connection with Landau-Zener-Stiickelberg 
(LZS) interference. P 13 E] 

Recently, a series of great progresses 0[TOl[IIl[I3[l3l[Tl[T5l[T6l[I71[T8l[l9l|2^ have 
been implemented to study the coherent dynamics of qubits in the regime dominated 
by LZS interference. The strong high-frequency microwave generates controllable 
population transitions, allowing for fast and reliable control of quantum systems. The 
delicate LZS interference patterns also provide a way to calibrate the parameters of 
the qubit and its interaction with the environment. On the other hand, through 
continuous Landau-Zener (LZ) transitions, the low-frequency microwave can be used 
to achieve an effective temperature of qubit up to two orders of magnitude lower 
than the bath temperature through an analogue of optical sideband cooling. [2T1 [22] 
This microwave cooling provides an effective means for preparing qubit state with 
high fidelity and suppressing decoherence. Therefore, understanding dynamics of 
LZS interference in strongly driven superconducting qubits has both fundamental and 
practical significances. A number of theoretical investigations [6l [HI [231 IMl [251 ISSl [271 
[28l[29l[3ni[3ll[3a[33l[ai[35l[Ml[371[38l^ have been developed in driven quantum 
two-level systems (TLS). These works generally addressed the coherent phenomena in 
the high-frequency driving in superconducting qubits, as reviewed recently in Ref. 6. 
In the realistic situation, the influence of environment always needs to be considered. 
Compared with natural atoms, the solid state systems contain macroscopic degrees of 
freedom strongly coupled to environment bathes that decohere the quantum states to 
be controlled. To address the driven dissipative TLS, two typical methods usually taken 
into account are rotating-wave approximation (RWA) [311 HO] and perturbation- RWA. 
[T9l[27] 

RWA provides a simple physical picture and good approximation for the stationary 
population. However, it limits at the high-frequency region, i.e. u ^ AJ„ (A/w), 
where n is the closest integer to Sq/uj and Sq is the dc detuning. The other method, 
perturbation- RWA, applies the second order perturbation approximation based on RWA 
and can be used in the parameter region with the high frequency and strong decoherence, 
i.e. w > A,r2 and r2 ^ VFoi? where Wm is microwave driven Landau-Zener (MDLZ) 
transition rate from |0) and However in the low-frequency driving where microwave 
cooling takes effects, both methods would be inappropriate. 

In this article, we focus on analyzing the stationary state which is useful in 
spectroscopy experiments [9l [TOl [121 IIS] and microwave cooling [21] in the presence 
of decoherence. Going beyond the second order perturbation and focusing on the 
nonresonant terms which are dropped in RWA, we present the real-time population 
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evolution as well as unsymmetrical MDLZ transitions in a strongly driven TLS of 
the long-time dynamics. In the high-frequency case, i.e. u > A, the result of 
perturbation approach can fully cover those of RWA and perturbation- RWA. In addition, 
the perturbation approach can be applied to larger parameters space, at which other 
methods are inappropriate. 

Then we move on lower frequency region, i.e. a; < A. In this parameter 
region, microwave cooling begins to take effects, motivating us to give a thorough 
understanding of the dynamics. The steady dynamics shows a periodic time dependence 
of qubit population, which is from the nonresonant terms and cannot be explained by 
perturbation-RWA. Then we focus on the average stationary population. Although 
one can use the perturbation approach, it is cumbersome thus difficult to get the 
final solutions. Therefore, we use the average perturbation approximation (APA), 
which is based on perturbation approach through the average over one driving period. 
Interestingly, the analytical expression in APA has the same form with that of 
perturbation-RWA. If the frequency further decreases, APA cannot correctly give the 
qubit population. We thereby propose the noncoherence approximation (NCA). The 
incoherent evolution is divided into discrete incoherent time steps, from which we can 
obtain the analytical solution. It is found that when the frequency is much smaller than 
the decoherence rate, the time dependent and average population calculated by NCA 
agree quantitatively with those of the numerical simulation. At last, we extend the 
analytical results of APA and NCA to the multi-level system. We discuss microwave 
cooling [21] and obtain the minimal frequency of the effective microwave cooling. 

Through these approximations of perturbation approach, we present simple and 
comprehensive descriptions in the average stationary population in most of parameter 
space of frequency and decoherence rate. The agreement with numerical results indicates 
that the picture captures the underlying physics, thus giving a better understanding of 
the behavior of system dynamics. The models and methods we used here are also valid 
for other TLSs with similar structures. 

The article is organized as follows. In section. 2, we introduce the basic model 
and approximations. In section. 3, assuming high order perturbation in the long-time 
dynamics, we present the perturbation approach which considers the nonresonant terms. 
Then with this approach we discuss the high-frequency driving case. In section. 4, We 
discuss the average stationary population in the low frequency. Through APA, we obtain 
the average population in the steady case under the low-frequency driving. In section. 
5, for the extremely low frequency, we demonstrate that NCA can be used to get the 
analytical solution, which is useful to describe the system dynamics. In section. 6, we 
extend our model to multi-level flux qubit and discuss the frequency range of effective 
microwave cooling. 
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Figure 1. (Color online), (a). Schematic picture of the double-well potential of a 
superconducting flux qubit. States |1) is in the right well; |0) is in the left well. (b). 
Schematic energy diagram of the flux qubit. Red solid curve represents the microwave 
^rf sinujt. The vertical dashed line marks a particular static flux detuning 5^dc- 



2. Basic model 

For a superconducting flux qubit [HI |12] with three Josephson tunnel junctions, if 
the external flux bias $ ^ 0.5$0) where $o = h/2e is the flux quantum, a double well 
landscape of the potential energy parameterized by the dc flux detuning 6^dc = 0.5$o 
exists as shown in Fig. 1. In the millikelvin temperature ~ 10 mK, a series of diabatic 
states are localized in different wells. States |0) and |1) are the lowest levels in the left 
and right wells respectively. They are characterized by different flowing directions of 
the persistent currents, namely clockwise and anticlockwise. The states in each well 
are separated by large energy spacing on the order of 20 GHz, [15] while the coupling 
strength between two wells is small compared with the energy spacing. As the qubit is 
driven with a microwave $ac = ^rf sin ut, where $r/ is the microwave amplitude, the 
time dependent flux detuning is 

= 5$^^ + = 5$^^ + sinut. 

Then the time dependent energy detuning of states |0) and |1) can be described as 

e{t) = + Asinut, 

where Eq = (|?7io| + |mi|)5$rfc is the dc energy detuning of states |0) and |1), A = 
(|mo| + |mi is the microwave energy amplitude, and rrii = dEi{^) / is the diabatic 
energy- level slope of state \i) in units of frequency per flux. For the initial time t = 0, 
the qubit is at the detuning e{0) = Eq- 

We start from a periodically driven TLS coupled to an environment bath. Such a 
process can be described by Hamiltonian (we set h = ks = I) [261 EZJ [39l HS} |44] 

A e(t) , , 

a = - -Y(r, + Hint + Hb, (1) 

where and are Pauli matrices, = |1) (1| — 10)(0| and = |1) (0| — 10)(1|, Hint is the 
system-environment interaction Hamiltonian, and Hb is the environment Hamiltonian. 
We assume the qubit is embedded in the bath in equilibrium at temperature T with 
the density matrix of environment pB = e~^^/^. The system-environment interaction 
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Figure 2. (Color online). Population in state |1) as a function of time, assuming 
that the system is initially in state In both figures, eo/27r = 90 MHz, A/27r = 90 
MHz, and A/2Tr ~ 5 GHz. The microwave frequency is (a) w/27r = 90 MHz and (b) 
ljj/2tt = 10 MHz, respectively . The blue solid line shows the numerical stimulations 
and the black dotted line shows the RWA results. 



Hamiltonian is written as 

H^ni = -Q^./2, (2) 

where Q is an operator acting on the environment. 

Before investigating the system dynamics in the presence of the environment bath, 
it is worth recalhng main features of its free time evolution which is described by 
Hamiltonian Eq. (1) without Hint and Hb- By a gauge transformation, the Hamiltonian 
is brought to the form 

H, = -^ef^'^\l){Q\ + H.c. (3) 

where (j){t) = —i j^e{T)dT. Using Bessel functions, we have e"'^*-*^ = 

With the relation e^^/^ = Jn{A/u)e'^"'"^'^, we consider parameters close to the n- 
photon resonance, where nu = Eq, and rely on RWA to neglect the fast oscillating terms 
with k n for the high frequency u. Therefore, when the system is initially in state |1) 
we have the population in state |0) [6l [9], [3l] 
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This approximation neglects all the nonresonant terms which oscillate very fast in 
the time scale of system dynamics, so the high frequency is needed to keep the validity 
of RWA, i.e. u > A ^ AJ„ (A/u). [HI EH [31] We compare the real-time numerical 
solution with the analytical results in Eq. (5), shown in Fig. 2. RWA has the high 
accuracy under the high-frequency driving. However, when the frequency is low, i.e. 
a; < A, the difference between numerical and RWA results are not negligible. It should 
be noted that in Fig. 2 (a) there exist two kinds of oscillations. One agrees with Eq. (5) 
while another as shown in the panel in Fig. 2 (a) is from the nonresonant terms in Eq. 
(4). When the frequency is small just as shown in Fig. 2 (b), the effect of nonresonant 
terms could not be neglected. In next sections, we will further show the influence of 
these nonresonant terms on system dynamics. 

In the realistic solid-state system, the decoherence which is always strong needs to 
be taken into account. In this case, the Bloch equations are usually used to provide the 
numerical stimulation, which read 

p = -z[H,,p]+T[p], (6) 

here T[p] is the dissipative term B6] 

■pr 1 — rioPii + ToiPoo — r2Poi 

= -p 

—J- 2P10 J- loPii — J- oiPoo 

where r2 = I/T2 is the decoherence rate, T2 is the decoherence time, Fjj is the interwell 
relaxation rate from state \i) to Fio = Fqic"^'^/"^, 1/Ti = maxjFio, Fqi}, and Ti is 
the interwell relaxation time. After the long-time system dynamical evolution, whatever 
the initial state is, the temporal oscillation in Fig. 2 is eliminated by the decoherence 
and the qubit population will converge to a stationary solution. We make RWA with 
the equation Eq. (6), neglect the fast oscillating terms, and obtain the probability in 
state |0) [6] 



(7) 



r i 10 

r^+(eo-nc^)2 + Toi + Tio 



To ensure Eq. (8) is valid, we need not only the high frequency, i.e. a; > A ^ AJ„ (^), 
but also the high coherence, i.e. u ^ If the decoherence is strong a; ~ the 
population poo niay even be larger than 1 and this stationary solution is obviously not 
reasonable. 

In order to describe the dynamics in the strong decoherence, Berns et al. [19] 
proposed the perturbation-rotating-wave approximation. Under the Gaussian white 
low-frequency noise, based on the second order expansion perturbation, we have MDLZ 
transition rate from state |0) to |1) 

_A'^ T,Ji(Alu) 
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Figure 3. (Color online). Stationary population in state |0) as a function of dc energy 
detuning in the presence of decoherence. In all figures, uj/2tt = 90 MHz, A/2tt = 90 
MHz, the interwell relaxation time Ti = l/(27r x 0.00005) ns, and the temperature 50 
mK. The driving amplitude is given by A/27r — 5 GHz in (a)-(c), and 0.01 GHz in 
(d)-(f). The decoherence rate is given by T2/2t: = 3 MHz in (a) and (d), 110 MHz 
in (b) and (e), and 1050 MHz in (c) and (f). The red solid line shows the analytical 
result from Eq. (11). The open circles give the numerical results. The inset shows 
clearly the accordance between analytical and numerical results. 



Correspondingly, Wiq^Eq) = Woi(— ^o) = Woi{eo). Then for the long-time dynamics, we 
employ the rate equations, in which the qubit population in state |0) and |1), i.e. poo 
and pii, follow 

Poo = -{Woi + ToOpoo + {Wio + rio)pii 

l = Poo + Pii- (10) 

Since MDLZ transition rate is constant in the stationary case, the equations would lead 
to the constant population, i.e. poo = Pii = 0, and the population in state |0) becomes 

AVg(A/^)r2 -p 

n 2^ ri+{eo-nwr + 01 + i 10 
n 

In order to ensure the validity of this approach. p!9| [27] the high frequency and strong 
decoherence, i.e. a; > A,r2 and r2 ^ Woi, are necessary, which greatly restrict the 
application of this approximation. In next sections, we will introduce approximations 
to extend the quantitative analysis to larger parameters space. 
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3. Perturbation approach 

In order to get analytical results in larger parameters space, we go back to Hamiltonian 
(1) and take the perturbation approach. First of all, we follow the process in Ref. 27 and 
44, and the complete time evolution operator in the laboratory picture can be expressed 
as 

U{h,t2) = e-*^«*^f/o(tl)f/l(^l)^/y(^l,^2)f/^(^2)f/o^(^2)e*''"*^ (12) 

where Uo{t) = T expff c/TC(r)(T,], Ui{t) = T exp[| J^* rfrg(r)a,], Q{t) = 
e'H^'Qe-^'^^Wviti.h) = rexp[-z/;;rfry(r)], V{t) = -| ^ J„(A/a;) exp[-z(£ot - 

n 

nut) + 17171/2 — iA/u\U\{t)\l){Q\Ui{t) + H.c, and T is the time ordering operator. Let 
PsB express the total density matrix of the system plus environment bath, which can be 
described as psBif) = U(t, 0)psb{0)W (t, 0). Then the system reduced density matrix is 
defined as p{t) = TiBipsBit)}, where Tr^l...} denotes the trace over the environmental 
degrees of freedom. We assume at the initial time the density matrix is separable, i.e. 
P5b(0) = p(0) ® Pb, and the system is initialized in state |0), i.e. p(0) = |0)(0|. 
We perform a perturbation expansion of f/y in V{t) , which reads 



f/y(t,0)^l-z / dTV{T)- I drVir) i dr'Vir') 
Jo Jo Jo 

+1 f drV{r) r dr'V{r') r dr"V{r") + ... (13) 
Jo Jo Jo 

For the strong decoherence, i.e. r2 ^ AJ„ (^) ~ Wqi, we can consider the 
expansion to the second order in a time interval l/r2 ^ t < l/[AJn{A/uj)]. |271 HI] 
However, in the high coherence region where r2 < AJn{A/uj), instead of considering 
the time interval t < l/[AJn{A/uj)], we need take account of a time interval t much 
larger than the decoherence time. Hence this time interval t is also much larger than 
l/[AJn{A/uj)]. Moreover, for the low frequency AJn{A/u) < w < A, in order to 
calculate MDLZ transition in the stationary case, even if l/[AJn{A/u))] is much longer 
than the decoherence time, the time interval t < 1/[A J„(A/a;)] is still not sufficiently 
long relative to the driving period. As a result we consider t ^ l/uj, l/[AJn{A/uj)], I/T2. 
It is clear that the second order perturbation expansion would be insufficient. [271 Si] 
Considering high order expansion terms, we obtain the population in state |1) (see 
appendix) 

Pii(^)~ r f drdr'TrB{{l\V{r)psB{0)V\r')\l)} 
Jo Jo 

-2/ / dhdt2TrB{{l\[ dtsdUV{U)psB{0)V\h)]V\h)V\h)\l)} 
Jo Jo Jo Jo 

-Tr-Bi dhdU{l\V{t,)\2 / dt2dhpsB{m\h)V\h)]V\h)\l)} + ... (14) 

Jo Jo Jo Jo 
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For the short time interval, i.e. t < l/[AJn{A/u)], we use the second order 
perturbation expansion and can keep only the first term in Eq. (14) expressing the 
population in state |1) transferred from |0), just as discussed in Ref. 27. As the time 
goes on, the population is not completely in state |0). The first term does not consider 
the effect of the reduced population in state |0) such that it cannot characterize the 
population transferred from state |0) to |1) any more. Higher order terms in Uv{t,0) 
need to be considered, which result into the second as well as third term in Eq. (14). 
In the second term, Jq Jq dt3dt4V{ti)psBiO)V'^ (t^) describes the population transferred 
from state |0) to |1) in the time interval t. Then the second term characterizes the 
transition of population from state |1) to |0) based on this part of population in 
|1). Similarly, 2 J^^ c?^2C^^3Ps_B(0)V^^(t3)V^^(t2) in the third term describes the reduced 
population in state |0) in the time interval t. The third term describes the decrease of 
population transferred from state |0) to |1) because of the reduced population in state 
|0). When the time interval t becomes further larger, more terms will appear in Eq. 
(14) as a ripple effect and we can also extend the integration limit t to oo. 

Then we only discuss the first term in Eq. (14) as an example, and other terms can 
be treated with a similar procedure. We rewrite the first term of Eq. (14) to 



y{t) = ^ / dh [ dhK{Ul{h)U+{h)UX{t2)U.{t2)) 

4 Jo Jo 
A A 

_ \ ^ J ^ ^ J- ^ '^^i(eo—rnijj)ti—i{eo—nuj)t2+imT/2-innT/2 

UJ UJ 



(15) 



where U±{t) = Texp[±(i/2) J^Q{T)dT]. For the low frequency, i.e. AJn{A/uj) < w < 
A, the terms in K with n ^ m would oscillate not fast enough, leading the invalidity of 
RWA. Assuming the white Gaussian low-frequency noise, we rewrite Eq. (15) to 

y{t) = — dt^ dhj^e^ 
^ Jo Jo ™ 



i{£0-'muj){ti-t2)-T2\ti-t2\ 



UU-) + E Jn+m{-)M-)e-'^^-/'}. (16) 
UJ UJ UJ 

For the long-time dynamics, we extend the integration limits of Eq. (16) to oo, and 
have [27J 

^i{£0—muj)T+inojt2-\-i'mT/2—r2 \t\ 



A"^ d f f A A 

fit) = — } Urn — dti dt2Jn+m{ — )Jm{ — 
4 Z-j t-^oo dt Jo Jo UJ u 



)e 

'--i>00 dt In In ^ Uj' ^ UJ' 

n^O,m 

where r = ti — ^2- The first term in Eq. (17) is the familiar result in perturbation-RWA. 
For convenience we change the form of f{t) to 

4 UJ UJ 

njtO,m 
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lim — f dr' f dre'^eo^"'-^)r+inuir'-r/2)+inn/2~r2\r\^ (^g^, 
t^oo dt Jo J_l 

where r' = {ti + t2)/2, and t = min{2r', 2t — 2r'}. Therefore we have the leading term 
given by 

f(f^^^' V 7 (^)T r2Cosn(a;t + 7r/2) 

Then we consider other terms in Eq. (14) and follow the above process, obtaining 

Pii(^) ~ y{t, Eq) - y{t, eo)piiit) - y{t, -eo)pii(i) 

= Woi{t)poo{t) - Wio{t)pii{t), (20) 

where Wqi (Wiq) is MDLZ transition rate from state |0) to |1) (|1) to |0)). We emphasize 
that in the derivation of Eq. (20), all the expansion terms in Eq. (14), not only the 
notable terms, have the contribution. Eq. (20) gives the evolution of the diagonal 
part of the system reduced density matrix under the "Markovian" approximation, (see 
appendix) It should be noted that we take the relaxation phenomenonally and hence 
it does not emerge in Hamiltonian and rate equations, while in the calculation of 
population we will add it. Furthermore, MDLZ transition rates obtained are the 
combination of Lorentzian line shapes and time dependent, 

2 ^-^ u! u (er, — muY + Ti 

n,m ^ ' ^ 

2 u uj [eQ — mujY + Ti 

n,m ^ ' ^ 

where we use the relation Jm{x) = (—1)"^ J^rnix) . Different from Eq. (9), MDLZ 
transition rates in Eq. (21) are unsymmetrical. Comparing Eq. (21) and Eq. (9), 
we find that their difference is the additional time-dependent oscillating terms f{t), 
which comes from the nonresonant terms. The nonresonant terms also correspond to 
the weak oscillation in Fig. 2 (b), due to the coherence evolution. It is clear that the 
coherence in the nonresonant terms results into unsymmetrical MDLZ transitions, i.e. 
Woi{t) 7^ Wio{t). In addition, although the coherence is preserved in the nonresonant 
terms, nonresonant terms with different oscillation frequencies counteract the resonant 
terms which would become less obvious. Therefore, not only the decoherence but also 
the low-frequency driving would eliminate the multiphoton resonant fringes. 
Then we have the population in state |0) 

Poo = e/o'^*'[-^oi(0-w^io(*')-ro,-rK,]| f dt'[W^,(t') + Fio] 



Similarly, by following the above procedure, we can extend perturbation approach 
to higher levels and different noise sources, which were previously explored under high- 
frequency driving. [261 123 Especially, MDLZ transition between higher levels can be 
used to analyze microwave cooling with the maximal efficiency. [22] 
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Under the large amplitude driving, the coherence of LZS interference comes from 
two aspects. The first one is traditionally from the prepared initial states. The off 
diagonal elements of the density matrix describe this part of coherence, which is usually 
needed to be preserved for the application of quantum computation. The second 
one is the coherence induced by the crossover which acts as a beam splitter. [9l [16] 
The microwave drives the state through the crossover, splitting the state into the 
superposition of states |0) and and the system gains the coherence again. For a 
stationary case, the first kind has disappeared fully while the second one is created 
incessantly. However, it should also be noted that usually the coherence induced by the 
crossover only exists between two consecutive passing through the crossover within one 
period. There is no coherence for the processes in arbitrary periods. 

To make this method reasonable, as discussed in appendix, the coherent evolution 
produced by the crossover cannot produce remarkable population change in the 
stationary case, which means AJ„(A/a;) has little effect in one driving period. If the 
tunnel coupling strength A is large and the decoherence is weak, the coherence will 
persist for a long time, producing remarkable coherent population change. Therefore, the 
validity of this method depends on two factors. The first one is the high frequency, i.e. u 
^ AJn{A/uj), which makes coherent dynamics not have enough time to cause population 
change. The second one is the strong decoherence, which makes the coherence disappear 
quickly and produce little effect on the population, i.e. r2 > AJn{A/u). If one of factors 
is satisfied, the method would be available. 

For the high frequency, i.e. w > A AJn{A/u), the time-dependent terms in 
Eq. (21) oscillate fast and can be neglected just like that in RWA. Therefore MDLZ 
transition Eq. (21) returns to symmetrical perturbation-RWA result Eq. (9) and is 
constant. In addition, as shown in the derivation, Eq. (9) can be further extended to the 
incoherent region u and the high coherence r2 < AJn{A/uj). Furthermore, in the 
condition, i.e. u S> T2, AJn{A/u), the stationary population Eq. (11) in perturbation- 
RWA is equivalent to Eq. (8) in RWA. To validate the perturbation-RWA result in 
the two regions, we have numerically solved the Bloch equations to compare with the 
analytical results obtained by Eq. (9) on the steady occupation probability in state |0), 
as illustrated in Fig. 3. Here we consider the frequency u/27r = 90 MHz with different 
decoherence rates. The analytical results agree well with the numerical simulations. 

For the low frequency, i.e. u < A, the nonresonant terms will be comparable with 
the oscillating terms and cannot be neglected. In this case, perturbation-RWA result is 
not available and Eq. (21) must be needed. As discussed above, if the decoherence rate is 
small, e.t. r2 < AJn{A/u) ~ Wqi, the strong coherent population evolution may break 
perturbation approach. Therefore in this region, for smaller frequency, perturbation 
approach can be used in stronger decoherence. When the decoherence is enough strong, 
i.e. r2 > AJn{A/uj), the coherence generated by the crossover would have no effects.. 

When the state is driven with the low frequency through the crossover A, the 
small sweeping velocity makes the effect of nonresonant terms released. Nonresonant 
terms generate the periodic change of MDLZ transition rate with the same period of 
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Figure 4. (Color online). Population in state |0) as a function of time in the presence 
of dccoherence, assuming that the system is initially in state |0). In all figures, 
£o/27r 4950 MHz. The amplitude is given by Ajl-K = 5 GHz in (a)-(c), 0.01 GHz 
in (d). The microwave frequency is given by w/27r = 10 MHz in (a), 90 MHz in (c), 
1 MHz in (b) and (d). The decoherence rate is T2/2tt = 110 MHz. Other parameters 
are identical with those in Fig. 3. The inset shows clearly the steady periodic change 
of population. 




the driving, which makes the interwell relaxation participate in the dynamics and cause 
the change of the qubit population, i.e. poo = Pii 7^ 0, in the stationary case. Then the 
population would also oscillate with the frequency u. As shown in Fig. 4, in the steady 
case with large amplitude A > eo, for the driving frequencies 10 MHz and 1 MHz, the 
population in state |0) exhibits a periodic change [Fig. 4 (a) and (b)] while for the 
frequency 90 MHz the population is always constant [Fig. 4 (c)]. On the other hand, at 
the small amplitude, the curve becomes smooth again [Fig. 4 (d)]. We can approximate 
AJn{A/u) = {A/n\){A/2u)"'. Small amplitude makes AJn{A/u) much smaller than u 
such that the high-frequency case Eq. (9) can be used. This case can also be understood 
simply: when the amphtude is small as illustrated, the state does not pass through the 
crossover and LZ transition would not occur, resulting little change of population in the 
steady case. 

4. Average perturbation approximation 

As discussed in Ref. 21 and 22, the low frequency u, i.e. AJn{A/u) < u < A, can 
produce the optimal microwave cooling which is used to realize the active cooling of 
qubit and other artificial devises through MDLZ transitions. Although the low frequency 
makes the qubit population time-dependent in the stationary case, we can consider the 
average population over time. For the high frequency, the average population is just 
the steady population described by Eq. (11) while for the low frequency, the average 
population can be obtained by perturbation approach. However, it is complicate and 
not convenient to calculate the average population with Eq. (21). Therefore, we will 
employ average perturbation approximation (APA) based on perturbation approach in 
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DC Energy Detuning 8q (GHz) 



Figure 5. (Color online). Average stationary population in state |0) as a function 
of do energy detuning in the presence of decoherence. In all figures, a;/27r = 1 MHz, 
and eo/27r = 4950 MHz. The driving amplitude is given by A/2Tr = 5 GHz in (a)-(c), 
= 0.01 GHz in (d)-(f). The decoherence rate is given by T2/2tt = 40 MHz in (a) and 
(d), 110 MHz in (b) and (e), and 1050 MHz in (c) and (f). Other parameters are 
identical with those in Fig. 3. The red solid line is the analytical result from Eq. (26). 
The open circles give the numerical results. 




DC Energy Detuning Sq (GHz) 

Figure 6. (Color online). Average stationary population in state |0) as a function 
of dc energy detuning in the presence of decoherence. In all figures, uj/2it ~ 10 MHz, 
and eo/27r ~ 4950 MHz. The driving amplitude is given by A/2Tr = 5 GHz in (a)-(d), 
= 0.01 GHz in (e)-(h). The decoherence rate is given by r2/27r = 40 MHz in (a) and 
(e), 110 MHz in (b) and (f), 1050 MHz in (c) and (g), 3 MHz in (d) and (h). Other 
parameters are identical with those in Fig. 3. The red solid line is the analytical result 
from Eq. (26). The open circles give the numerical results. 
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this section. 

Then, we go back to the expression of f{t), Eq. (19). In the stationary case, 
considering the time in one period, i.e. kT < t < kT + T, for the term M{n, m, t) in 
fit), we have 

M n, m, t) = Jn+m{-)Jm{-)^ S^^' ^ ^ 0' 23 

where M{n,m, kT + t') = M{n,m, kT + T/2n + t') with an even number n and 
M(n, m, kT + t') = M{n, m, kT + T — t') with an odd number n. Therefore, in one 



l-kT+T 
kT ' 

average MDLZ transition rates by averaging Eq. (21) over one period in the stationary 



period of time, the above relation leads to lim J^^ dt' fit') = 0, and we obtain the 



case 



1 pkT+T 

{W,o) = {Woi) ^ - lim / dt'Woiit') 

_A^^ T,4iAM 

2 2^{e,-nuj)^ + rl- ^ ' 

n 

It should be mentioned that although Eq. (24) coincides with the transition rate 
Eq. (9) in perturbation-RWA, the latter is the stationary transition rate under the 
high-frequency driving while the other is the average transition rate through averaging 
over the driving period under the low-frequency driving. Having obtained the transition 
rate, we can describe the population of system dynamics using the rate equations, where 
the average population (pjj) obey 

(poo) = -((Woi) + roi)(poo) + ((W^io) + rio)(pii) 

1 = (poo) + (Pii) (25) 

After the long-time system evolution, we have (poo) = (pii) = 0, and the average 
qubit population in state |0) can be easily solved from Eq. (25) 

(poo) = ((W^oi) + rio)/(2(W^oi) + Tio + Toi). (26) 

The rate equations used here are different from Eq. (20). Under the low- 
frequency driving, Eq. (20) describes the real-time dynamics with the periodic change 
of population. However, Eq. (26) neglects the population change within one period 
as shown in Figs. 4 (a) and (b) and focuses on the scale of one period, in which the 
average population would be constant. If APA is valid, the change of qubit population 
must be slow on the scale of the driving period, i.e. uj > (Wqi). If u is larger 

t' 



than {Wqi), the population change in one period is so small that j^,^^ WQip^odt/T 



{poo) J^,~^^ Woidt/T = (poo)(W^oi) and the rate equations can be changed to the form of 
Eq. (25). 

In order to show the region of validity for APA, we compare the numerical solutions 
with analytical results. The average transition rate (Wqi) is probably less than 5 x 27r 
MHz and generally the same order with AJn{A/u). For several decoherence rates, we 
select the frequency 1 MHz and obtain the average steady population, shown in Fig. 
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5. The deviation between the numerical and analytical results cannot be neglected 
and APA is not proper. Fig. 4 (b) shows the real-time population in state |0) with 
the numerical stimulation under the frequency 1 MHz. Comparing with the result of 
frequency 10 MHz in Fig. 4 (a), due to u < (Wqi) the population change in Fig. 4 (b) 
is large, demonstrating the invalidity of APA. (Fig. 5) 

Moreover, as shown in Fig. 5 (a-c), at the frequency 1 MHz, when the amplitude 
is less than the dc detuning, which means the state does not reach the crossover, APA 
exhibits the good agreement with the numerical results. This is more clearly seen in 
Fig. 5 (d-f) with a small amplitude A/2n = 10 MHz. Since the population change in 
the steady case is little [Fig. 4 (d)], APA is still available. In Fig. 6, for the frequency 
10 MHz, the analytical results show very good agreement for several decoherence rates, 
changing from the rate less than AJn{A/u) to the one much larger than A. When the 
amplitude is larger than the dc detuning, the population is near 0.5, which also provides 
a way to determine the location of crossover with the low frequency. [Fig. 6 (a-d)] It 
should be noted that with the decoherence T2/27T = 3 MHz and the frequency 10 MHz, 
the agreement is not very good at the detuning Eq ~ A. The high coherence makes 
the average transition rate large. As discussed in Ref. 22, at the detuning Eq ~ A, 
the transition rate reaches maximum. Therefore the frequency would be less than the 
average transition rate, i.e. u < (Wqi), and results into the disagreement. 



5. Noncoherence approximation 

For the extremely low frequency, i.e. u < AJn{A/u)) ~ (Wqi), APA becomes insufficient 
just as shown in Fig. 5. Therefore in order to analyze this region conveniently, here we 
develop an analytical method based on noncoherence approximation (NCA), which will 
give an intuitive and clear physical picture. 

When the driving frequency is high and coherent, LZS interference and multiphoton 
effect are clear and Bessel functions can give a convenient picture. However, the decrease 
of frequency makes LZS interference gradually disappear. Therefore the description with 
Bessel functions would not be proper for the extremely low frequency and we return back 
to the sinusoidal description. 

We go back to Eq. (16) and have 

f.r'+r/2 



d /■* /■* rr'+T/2 

y{t) = —— dr' dTexp{i dt' {eq + Asinut') - T2\t\} 

4 dt Jq J_i Jt'-t/2 

A^ /•* /•*+^ 
= — / dTexp{i / dt'{Eo + Asinujt') -T2\t\} 
4 Jo Jt 

A2 /■* 

+— / drexpii dt'{Eo + A sin cot') - T2\t\}. (27) 

4 J-t Jt-T 

Because of the strong decoherence, the extreme low frequency usually is much less 
than the decoherence rate. The decoherence makes the qubit lose coherence quickly, even 
though the state is split by the crossover and obtains the coherence again in each period. 
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Hence we can approximately replace smut' to sin cut. Then we follow the procedure in 
Sec. 3, and have MDLZ transition rates 

W^t) = WoM = (28) 

Eq. (28) is actually the result of Eq. (21) in NCA limit. Interestingly, MDLZ 
transition rates become symmetrical. The expression can also be understood as 
macroscopic resonant tunneling (MRT) ^31 H5l W7\ model under the sinusoidal driving 
with the extremely low frequency. When the state is driven slowly to pass through each 
detuning, we divide the evolution into discrete time steps within which the coherence 
lost totally at the corresponding detunings. [H] Furthermore there is no coherence to 
connect time steps and the phase accumulation [9] has not existed. Since each time 
step is much larger than the decoherence time, we can approximately consider that 
MRT occurs at the detuning with the tunneling transition rate Eq. (28). Hence, the 
dynamics is changed from LZS interference to the periodically driven MRT and we have 
the rate equations 

Poo = -[Woiit) + roi]poo + [Wio{t) + Tiolpii 

l = Poo + Pii, (29) 

where Wqi is defined by Eq. (28) and roi(rio) has the same definition with that in 
Eq. (7). In the evolution equations Eq. (29), the periodic change makes the previous 
treatment poo = pn = incorrect. When the initial state is |0), we have the population 
in state |0) 

Poo = e/o'rf*'[-2H^oi(t')-roi-r,o]| f dt'[Wi,{t') + T,o]e^">'"^-^'^^^^'"^-^^^-^^-^ + 1}. (30) 

Jo 

The limit of perturbation approach also has function on NCA, which requires the strong 
decoherence, r2 > AJ„(A/u;) ~ (Wm), where {Wqi) is defined by Eq. (24). Moreover 
NCA needs the extremely low frequency compared with the strong decoherence. 
Therefore, in this region, this method is reasonable with the condition u < (Wqi) < 
In the frequencies where the average perturbation theory fails, we use Eq. (30) to 
reconsider the average population in the steady situation as shown in Fig. 7. At the 
decoherence rate r2/27r = 3 MHz, the agreement is not good, because the decoherence 
rate is comparable with the driving frequency. However, for other decoherence rates, the 
numerical stimulations agree well with the analytical results. Furthermore, we consider 
the real-time evolution and compare it with the numerical stimulation (Fig. 8). As the 
frequency is comparable or larger than the decoherence rate, NCA will be insufficient. 
[Fig. 8 (b)] 

From perturbation-RWA to NCA, we can find the transition from quantum to 
classic effect of the microwave with the frequency decreasing. For the high-frequency 
driving, being in a coherent state and with high photon numbers, the microwave behaves 
as the quantized electromagnetic field exchanging photons with the qubit. At the same 
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Figure 7. (Color online). Average stationary population in state |0) as a function of 
dc energy detuning in the presence of decoherence. In all figures, uj/2-k = 1 MHz, and 
A/27r = 5 GHz. The decoherence rate is T2/2'k = 3 MHz in (a), 40 MHz in (b), 110 
MHz in (c), and 1050 MHz in (d). Other parameters are identical with those in Fig. 
3. The red solid line shows the analytical result in Eq. (30) through noncoherence 
approximation. The open circles give the numerical results. 
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Figure 8. (Color online). Population in state |0) as a function of time in the presence 
of decoherence, assuming that the system is initially in the state |0). In both figures, 
Lo/i-K = 1 MHz, £o/27r = 4125 MHz, and A/2'n = 5 GHz. The decoherence rate is given 
by r2/27r = 1050 MHz in (a), and 3 MHz in (b). Other parameters are identical with 
those in Fig. 3. The red line shows the results from the noncoherence approximation. 
The open circles and dotted line give the numerical results. 
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Figure 9. (Color online), (a) Population in state |0) of TLS described by Fig. 1 
as functions of microwave frequency and decoherence rate. Eo/Stt ~ 4950 MHz, and 
Ajl-K = 5 GHz. Other parameters are the same with those in Fig. 3. The red 
surface with dashed line is obtained from the average perturbation approximation. 
The yellow surface with solid line is resulted from the noncoherence approximation. 
The blue solid circles are from the data of numerical stimulation, (b). The detail of 
(a) at the decoherence rate V^j'^.-n = 110 MHz. (c). Regions of validity for different 
approximations. The axes are the driving frequency and decoherence rate, respectively. 
Perturbation-RWA region is described by the condition a; > A and it is marked by 
the red circles. The average perturbation approximation region is described by the 
condition (W^io) < w < A and it is represented by the blue squares. The noncoherence 
approximation region is shown by the condition uj < (M^io) < r2 and it is marked by 
the black ellipses. 
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time, the large microwave amplitude forms LZS interferometry through the crossover, 
modulating the intensity of the n-photon resonance. In this case MDLZ transition Eq. 
(9) is based on perturbation-RWA and the nonresonant terms in the transition rate 
can be neglected. With the frequency further decreasing, the decoherence becomes 
stronger relative to the driving frequency and the nonresonant terms are comparable 
with the resonant ones, which makes the photons resonances as well as the modulation 
from LZS interferometry vanish gradually. Then we have to use Eq. (21) to describe 
MDLZ transition. For the low frequency, in the scale of one driving period, the resonant 
information can still be notable and we can use APA to obtain the average population. 
For the extremely low frequency compared with the decoherence, the microwave can 
finally be treated as a classic field and the dynamic process is changed to the periodically 
driven MRT as shown in Eq. (30). 

To sum up, we plot the average stationary population obtained from analytical 
calculation and numerical stimulation as shown in Fig. 9 (a). The different 
approximations of perturbation approach cover most of parameter space. Fig. 9 (b) 
is the detail of Fig. 9 (a) at the decoherence rate T 2/211 = 110 MHz. We show the 
applicable regions for these approximations in the average stationary population in Fig. 
9 (c). The circles mark where perturbation-RWA holds. In the red shaded region I, 
a high-frequency coherent field is used to realize multi-photon resonance. [9l [19] With 
the frequency further decreasing, we enter (Woi) < a; < A. In the stationary state 
the qubit population exhibits the periodic oscillation which cannot be explained using 
perturbation-RWA. Now APA is applicable, marked with blue squares in Fig. 9 (c). 
In the green shaded region II, the low-frequency driving is able to generate microwave 
cooling and the optimal microwave cooling was realized by Valenzuela et aL.[21] Then 
for the extremely low frequency u < (Wqi) < APA is also unsuccessful. NCA 
can describe the behavior, as marked by the black ellipses. In the blue shaded region 
III, which is not covered by any before mentioned approximation, the high coherence 
and low frequency make it possible and convenient to realize adiabatic LZ process or 
LZS interference by a single crossing [161 SHI ED] with large amplitude. It should be 
emphasized that the region in Fig. 9 (c) would change with the amplitude just as shown 
in Fig. 5. For the amplitude larger than the dc detuning, the small frequency locates 
in NCA region while for the amplitude less than the dc detuning, this frequency may 
move to APA region. 

6. Multi-level system 

Our approach gives a simple but effective physical picture, which enables us easily to 
extend the discussion to a multi-level flux qubit [Fig. 10 (a)]. [151 EI] In this section, 
based on the approximations we will discuss the lower limit of driving frequency in 
effective microwave driven cooling. 

In Ref. 21, it was found that a lower driving frequency realized better cooling effect 
and at 5 MHz the optimal cooling appeared. The driving microwave frequency used is 
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Figure 10. (Color online), (a). Schematic energy diagram of a multi- level flux qubit. 
Red solid curve represents the microwave ^ac- The dashed black line marks a particular 
static flux detuning S^dc- State |0) is in the left well; |1) and |2) are in the right well. 
The red solid path describes the cooling with the population in state |0) transferred to 
|1). The locations of the crossovers Aqi, and A20 are m<I>o, and 8.4 m$0 7 respectively. 
The diabatic energy-level slope |?no|(|wi|) =2t: x 1.44 GHz/m<i>o, and \m2\ =27r x 1.09 
GHz/m$o- Aoi/27r = 0.013 GHz, Ao2/27r = 0.09 GHz. The parameters used above are 
from the experiment. [15] In calculation, r2/27r = 0.06 GHz, r2i/27r = 0.1 GHz, [27] and 
the interwell relaxation time Ti = l/(27r x 0.00005) ns. (b). The yellow surface with 
solid line shows the population in state |0) versus microwave frequency and amplitude 
at dc flux detuning 6^dc = 0.05 m$o- The temperature is 50 mK. The dashed line 
indicates the population in equilibrium at this detuning. 



SO low that one does not need the high-frequency microwave generator and microwave 
hne to reahze active coohng. On the other hand, generally, if the microwave frequency 
is too low, there would not be the cooling effect on the qubit. Therefore, there exists a 
minimum frequency to produce the effective microwave cooling. Considering the three- 
level system consisting of an "internal" qubit (|1) and |0)) and an "internal" oscillator 
like state |2), as shown in Fig. 10 (a), we rewrite the rate equations Eq. (29) to 
three-level form 

Poo = -[W^Ol(t) + Toi + W^02(t)]p00 + [W^io(t) + Tiolpil + W2o(t)p22 

Pii = [W^oi(i) + ToiIpoo - [W"io(i) + Tiolpii + r2lP22 

1 = Poo + P1I +P22- (31) 

where is the intrawell relaxation from state |2) to |1). 
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For the low frequency as discussed in Sec. 4, APA can be applicable and we have 
W it)-W (t) - V (r2 + r.i/2)J^(A/a;) 

and, 

E^i^^!^. (33) 

where = {\mQ\ + \m2\){S^dc — ^2o), = (l^^ol + I'^iD'^^dc and $20 is the flux detuning 
where the crossover A20 is reached. It should be mentioned that Eqs. (32) and (33) 
actually express the average transition rate. Then for the extremely low frequency as 
discussed in Sec. 5, we have 

W,,it) = Wo.it) = ^ ^J't^':'l (34) 



2 e*{ty + {T2 + T2i/2y 



and, 



H'o.W = W.„(t) = ^^^^^ (35) 

where e{t) = (|mo| + |mi|)(5$(t), and e*{t) = {\mo\ + \m2\)[6^{t) - $20]- Here, we 
focus the extremely low frequency region. Since the analytical solution is too complex 
to extract a clear physical picture, we do not write it out explicitly here. Fig. 10 (b) 
shows the population in state |0) at the dc flux detuning 6^dc = 0.05 m$o as functions 
of frequency and amplitude. The thick dashed line marks the population at equilibrium 
with the temperature 50 mK. It can be seen that when the frequency is less than 0.1 
MHz, microwave cooling would be ineffective. 

In the extremely low frequency, the interwell relaxation changes the qubit 
population within one period as shown in Eq. (31). When the state is driven through 
the crossover A20, LZ transition would transfer the population from state |0) to |2), 
which then quickly relaxes into state |1). As the states is driven back to the initial 
detuning, the interwell relaxation warms the qubit again if the driving frequency is less 
than the repopulation rate of the warm process, i.e. Fio + Foi, which is approximately 
2Fio at the small dc detuning. The analysis agrees with the result in Fig. 10 (b). 



7. Conclusion 



In this article we have studied the population evolution of a strongly driven 
superconducting qubit of the long-time dynamics. First of all, we use perturbation 
approach, which considers high order perturbation expansion and nonresonant terms, 
to obtain the time dependent unsymmetrical MDLZ transitions. Moreover, we point 
out that this approach has a high validity in the case where the influence of coherence 
resulted from A is small with the high frequency or strong decoherence. Then for 
different frequencies, we apply several approximations to discuss the system dynamics. 

For the high frequency, i.e. a; > A, MDLZ transition rate in perturbation 
approach is constant and leads to perturbation- RWA result. In previous works [9], [19] 
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perturbation-RWA has been proven successful in the coherent region under the strong 
decoherence, i.e. u > A,T2 and r2 ^ Wqi. Besides this region, we show perturbation- 
RWA result can also be used in the incoherent case where the frequency is smaller than 
the decoherence rate, r2 > u, and even the high coherence region r2 < VFqi- 

With the frequency further decreasing to the region of the optimal microwave 
cooling, (Wio) < a; < A, since the coherence has effect on the nonresonant terms 
and nonresonant terms are comparable with the resonant terms, the qubit population 
exhibits periodic oscillation which cannot be addressed by perturbation-RWA. In this 
case, we study the average stationary population in the low-frequency driving. Then 
we employ APA by averaging MDLZ transition rate of perturbation approach with 
time. Although the analytical expression is the same with that of perturbation-RWA, 
they have different physical pictures. In perturbation-RWA, nonresonant terms have no 
effects and can be neglected, while in APA the nonresonant terms result oscillations and 
can be averaged over to obtain the average population. Therefore in APA, we do not 
take account into the real-time dynamics, but consider one period as the unit time. Then 
in the extremely low frequency, i.e. u < (Wiq) < where APA has failed, we propose 
NCA by dividing the evolution into discrete time steps during which the coherence 
is lost totally. We show the dynamic process equivalent to the periodically driven 
MRT, giving a convenient and intuitive picture. In addition, in the multi-level system, 
based on this approximation we examine the possibility of employing the extremely low 
frequency driving as a useful tool to realize microwave cooling and prove that there exists 
a lower limit of frequency for effective microwave cooling. The different approximations 
of perturbation approach cover most of parameter space when describing the average 
population. [Fig. 9 (c)] 

Our theory offers the effective analytical description of the driven dissipative TLS 
in the stationary case, including real-time and average situation. The good agreement 
of our theory with the numerical results shows the validity of the approximations used. 
By shifting these approximations, we also show the transition from quantum regime 
to classic regime, with dynamics changing from LZS interference to periodically driven 
MRT. This could enhance the understanding to the influence of decoherence on the 
qubit. Furthermore, our results can be applied to understand the behavior of other 
systems with similar structure. Especially, our theory in the low frequency can be used 
to analyze the optimal microwave cooling of qubit as well as other artificial device such 
as nano mechanical resonators. 
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Appendix A. Derivation of Eq. (14) 

In this appendix, we will accept tlie straiglitforward pliysical picture for intuitive 
understanding. Assuming that the system is initialized in state |0), first of all, we 
try to calculate the population in state |0), i.e. poo(t). In a short time interval, i.e. 
t < l/[AJn{A/uj)], considering the perturbation expansion of Uy in V{t) up to the 
second order term, we have 



Pooit) =TrBmpsB{tm} = TrBmUv{t,0)psB{0)Ul{t,0)\0)} 



t PT 



Jo 



TrBmpsB{0)m-TrB{ / drdr'(0|V(r)V(r')p5B(0)|0)} 



-TrB{ I r drdr'{0\psB{0)V\T')V\r)\0)}. (A.l) 
Jo Jo 

We substitute V(t) into Eq. (A.l) and have 
Jo Jo 

r drdr'K^iJJl {t')U+{t')uI{t)U.{t)) 
Jo 

= 1-2/" r dTdT'Ki{Ul{r)U+{r)Ul{T')U_{T')), (A.2) 
Jo Jo 

where 

A A 

\ ^ J (_\J ( \^i(eo—mijj)T—i(eo—nui)T'+inn/2—imTr/2 



and 



\ ^ J {^iS^J ^ ^ i(£o— mi.j)r+i(eo— na;)r'— m7r/2+)im,7r/2 



In the above derivation, we also obtain TrB{{^\ drdr'V {t)V {t')psb{^)\^)} = 
Trsm jlrodrdr'psBmKr'WKrm}- Moreover 2Tr^{(0| To drdr'V {r)V {r')psB{m)} 
just describes the reduced population in state |0) in the short time interval. 

Then we turn to high order expansion. It should be mentioned for simpleness in 
the below derivation we only consider the first four terms in Eq. (13) and higher order 
terms can be treated with the similar process. We have 



PnW =TrB{(l|f/y(t,0)p5s(0)f/i.(t,0)|l)} 
dTdT'TrB{{l\V{r)psB{mKr')\l)} 



t pt 



JO 

t ft pti l-t2 

dUdtidt2dhTrB{{l\V{U)psB{m\h)V\h)V\t{)\l)} 



Jo Jo Jo 



'0 ^0 JO JO 



dtidt2dhdt4TrB{{l\V{h)V{t2)V{h)psB{0)V\t4)\l)} + ... (A.3) 
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Considering that the last two terms in Eq. (A. 3) is equivalent, we can write Eq. (A. 3) 



as 



PllW= r f dTdT'TrB{{l\V{T)psBiO)V\T')\l)} 

Jo Jo 

-2 [ [ dUdhdt2dhTrB{{l\ViU)psBiO)V^{t3)V\t2)V\h)\l)} 
Jo Jo Jo Jo 



Jo JO 

' drdr'TrB{{l\V{T)psB{0)V\T')\l)} 

^0 

t fti l't2 pt2 



/*t /*^2 /*^2 

-2 / / / / dhdt2dtsdUTrB{{l\V{U)psB{0)V\ts)V^{h)V\U)\l)} 
Jo Jo Jo Jo 

-2 I ! dhdUdt2dhTrB{{l\ViU)psBiO)V^{t3)V\h)V\h)\l)} + ... (A.4) 

Jo Jt2 Jo Jo 

Then we find that the system population rate change is given by 

Pii(t) = f dT'TrB{{l\V{t)psBm^{T')\l)}+ f drTrB{{l\V{T)psB{0)V\t)\l)} 
Jo Jo 

-2 [ r [ ' dt2dt3dUTrB{{l\V{U)psB{0)V\ts)V\h)V\t)\l)} 
Jo Jo Jo 

-2 [ r f dt2dhdUTrB{{l\V{U)psB{m\h)V\t2)V\t)\l)} 

Jo Jo Jt2 

-2 [ r r dhdt2dhTrB{{l\V{t)psB{m\h)V\h)V\h)\l)} + ... (A.5) 
Jo Jo Jo 

Combining the equivalent terms in Eq. (A.5), we have 

Pii{t)=TrB{2 f dT{l\V{t)psB{m\r)\l)} 
Jo 

-TrB{2 [ dt2{l\[f' [ ' dt-sdUV{U)psB{0)V\h)]V\h)V\t)\l)} 
Jo Jo Jo 

-TrB{2 I dt^{l\V{t)[2 I ' dt2dhpsB{mKh)y\t2)]VKh)\l)} + ... (A.6) 
Jo Jo Jo 

where j^^ j'l^ dt3dtiV{ti)psB{^)V\tz) describes the population transferred from state 
|0) to |1) in the time interval ^2- Comparing with Eqs. (A.l) and (A. 2), we can 
find that the second term in Eq. (A.6) characterizes the population rate transferred 
from state |1) to |0) based on the population transferred from state |0). Similarly 
2 dt2dt3psB{^)y\tz)y\t2) describes the reduced population in state |0) in the 

time interval ti and the third term in Eq. (A.6) describes the reduced population 
rate from state |0) to |1) due to the decrease of population in state |0). Then we 
replace the integration limit ^2 in Eq- (A.6) to t, which is similar with the widely known 
Markovian approximation dealing with the population transition caused by the coupling 
to the bath [51]. Here, different from the Markovian approximation, the approximation 
is used in the population transition caused by the tunnel coupling A. To make this 
approximation valid, the system should be independent on the past, which means the 
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coherent population change resulted from A should be little. Integrating with Eq. (A. 6) 
we can obtain Eq. (14). 
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